Quantum Criticality in a Bosonic Josephson Junction 
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In this paper we consider a bosonic Josephson junction described by a two-mode Bose-Hubbard 
model, and we thoroughly analyze a quantum phase transition occurring in the system in the limit of 
infinite bosonic population. We discuss the relation between this quantum phase transition and the 
dynamical bifurcation occurring in the spectrum of the Discrete Self Trapping equations describing 
the system at the semiclassical level. In particular, we identify five regimes depending on the strength 
of the effective interaction among bosons, and study the finite-size effects arising from the finiteness 
of the bosonic population. We devote a special attention to the critical regime which reduces to the 
dynamical bifurcation point in the thermodynamic limit of infinite bosonic population. Specifically, 
we highlight an anomalous scaling in the population imbalance between the two wells of the trapping 
potential, as well as in two quantities borrowed from Quantum Information Theory, i.e. the entropy 
of entanglement and the ground-state fidelity. Our analysis is not limited to the zero temperature 
case, but considers thermal effects as well. 



INTRODUCTION 



Several years ago it was suggested that two noninter- 
acting condensates trapped in a double well trap at zero 
temperature would give rise to the bosonic analog of a 
Josephson junction Further studies of such a bosonic 
Josephson junction followed, where the effect of boson in- 
teraction was taken into account in a semiclassical frame- 
work based on the so-called discretized Gross-Pitaevskii 
equations '^-^ . Different interesting regimes and behav- 
iors were predicted for the system depending on the value 
of the effective interaction among the bosons in each of 
the two wells of the trapping potential. These include 
the so-called macroscopic self-trapping phenomenon and 
a dynamical bifurcation corresponding to a localization 
transition in the ground-state of the system Q. Just a 
few years after the above mentioned theoretical work, 
brilliant experiments were carried out with ultra-cold 
atomic gases where the predicted Josephson oscillations, 
macroscopic nonlinear self-trapping and dynamical bifur- 
cation were observed 041l- 

The quantum version of the discrete Gross-Pitaevskii 
equations — known as discrete self-trapping (DST) equa- 
tions [l^l among the nonlinear community — was also 
widely investigated [Tl|-[l3l. In particular, a significant 
amount of attention lj2j| was directed to the quan- 
tum counterpart of the abovementioned dynamical tran- 
sition from a localized to a delocalized regime, which, in 
this paper will be clearly recognized, as a quantum phase 
transition. 

In particular we carry out a detailed analysis of the 
ground state of the system which, at the quantum level, 
is described by a two-mode Bose-Hubbard model. We 
identify five regimes depending on the effective bosonic 
interaction. The central, critical regime corresponds to 
the dynamical bifurcation occurring in the solutions of 
the time-independent DST equations at the semiclassi- 



cal level. Actually, this regime shrinks with increasing 
bosonic population, and reduces to the semiclassical crit- 
ical point in the limit of infinite population which, in 
the present setting, has the role of a thermodynamic 
limit 2J]. Focusing on the "standard" order parame- 
ter, i.e. the population imbalance of the two wells, as 
well as on quantities borrowed from quantum informa- 
tion theory, such as the entropy of entanglement and 
the ground-state fidelity , we calculate the critical expo- 
nents characterizing the quantum phase transition in the 
thermodynamic limit. Moreover, we study the finite-size 
scaling effects, evidencing, within the critical regime, the 
anomalous scaling of the considered quantities as a func- 
tion of the total boson population. Our analysis is not 
restricted to the ground-state of the system, i.e. to the 
zcro-tcmpcrature case. We also consider finite temper- 
atures, and evidence a different scaling of the quantum 
and thermal fluctuations of the population imbalance in 
the critical regime. 

The plan of the paper is as follows. In Section |lT] we 
describe the quantum Hamiltonian employed in the de- 
scription of the two-mode system, and discuss its relation 
with similar models, as well as with the semiclassical DST 
equations; In Section lllll we illustrate how the low-lying 
spectrum of the quantum Hamiltonian can be equiva- 
lently investigated by means of a Schrodinger-like equa- 
tion for a flctitious particle on a one-dimensional domain, 
and observe that the semiclassical DST equations natu- 
rally emerge from a small-oscillation approach around the 
minima of the potential energy. Section |IV] discusses the 
so-called Gaussian regimes, where the Schrodinger-like 
equation can be safely approximated by a Schrodinger 
equation proper, featuring a harmonic potential part. 
The limits of validity of such an approximation are ana- 
lyzed in Section |Vl where further non-Gaussian regimes 
are discussed. Some more technical issues concerning the 
strong-coupling regimes are deferred to Appendix [Xj The 
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critical regime is detailedly analyzed in Sec. I VII whereas 
Sec. IVIII is devoted to a discussion of the behavior of 
two interesting indicators borrowed from Quantum In- 
formation Theory, the entropy of entanglement and the 
ground-state fidelity. The thermal fiuctuations of the or- 
der parameter are investigated in Sec. I Villi while Sec lIXI 
contains our conclusions and perspectives. 



is a su(2) coherent state, as discussed in Ref. [31 . The 
dynamics of the above variables is determined by the 
expectation value of the trial state on Eq. ([T]), which 
plays the role of a semiclassical Hamiltonian, 

n= mH\^) (^^z^ ~ ez- ^/T^^cosipj (6) 



II. THE MODEL 



The Hamiltonian of the system reads 

_u n V 



(1) 



where 



-ffint = ( a\a\aiai + a\a\a2a2 



= [ni(ni - 1) + n2{n2 - 1)] 
ffkin = {a\a2 + 4ai) 
-f^bias = {ni - n2) 



(2) 
(3) 
(4) 



and the lattice operators a 



and 



^. cLj.ua itj — ujjUij respec- 
tively create, destroy and count bosons at lattice site j. 
The operators in Eqs. ^ and ^ account for interac- 
tions among bosons at the same site and for their kinetic 
energy, respectively. The relative importance of these 
terms is controlled by the interaction strength U and the 
hopping amplitude 51 > 0. The operator ifbias allows for 
a local energy oflFset between the two sites, breaking the 
mirror symmetry of the system. 

Notice that Hamiltonian ^ commutes with the total 
number operator, iV = ni + 712, which therefore can be 
considered as a constant for all practical purposes. Owing 
to this fact, the inclusion of a further term of the form 
r^fdip = Tnin2, accounting for an interaction between 
bosons at different lattice produces a Hamiltonian of the 



same class as H 

Also notice that 
tween the spectral 
for attractive U < 
teractions. This 
serving that 



251 ■ Note indeed that H. 



dip 



a simple mapping exists be- 
features of Hamiltonian ([T]) 
and repulsive U > in- 
can be appreciated by ob- 

D [^Hiat + f i?bias — §ffkin] = 

- [-f ifint - f ^fbias -_|i?ki„] , where D = e*-(»i-"2) is 
a unitary transform, |2(f . 

Hamiltonian ([T]) can be studied by assuming that the 
state of the system is well approximated by a trial wave- 
function of the form 




|0) (5) 



where — 1 < ^ < 1 and ip are two real dynamical variables 
describing the population imbalance and the macroscopic 
phase difference between the two wells [1, [l^j . Eq. (O 



where 7 = and £ — fi are the effective param- 
eters of the model [Ij. Specifically, the dynamics is 
dictated by the so-called Discrete Self-Trapping equa- 
tions ensuing from Eq. and the corresponding 
fixed-point equations are known to exhibit bifurcations 
at |7|2/3 - 1 = |e|2/3 [i§ [2|. The manifestation of the 
semiclassical bifurcations in the features of the quantum 
system has been repeatedly reported in the literature. 
For instance the emergence of structures in the spectrum 
of the quantum problem ([T]) at the same energies as the 
semiclassical fixed- poin t solutions has been investigated 
in Refs. [H, M, 

In the following we will be mainly concerned with the 
mirror-symmetric case, e = 0, for which the bifurca- 
tion corresponds to a dynamical transition at 7 = — 1 
in the ground-state of Hamiltonian For 7 > — 1 

this is symmetric, z = 0, while for 7 < — 1, due to the 
nonlinear nature of the dynamical equations, the mirror 
symmetry of the Hamiltonian is spontaneously broken, 
z = iy'l -7-2. 

Although the ground-state of the quantum Hamilto- 
nian ([1]) is always symmetric, (ni) = (112), a crossover 
in its internal structure can be recognized in the vicin- 
ity of the mean-field critical value. Several indicators 
have been employed to highlight this crossover: energy 



fT4,l3lll32l | , flu ctuations in the population imbalance 
1 Il6l 19l. l2Cll. I23I . structure of the occupation number 



fidelity [31[, entanglement 
Hliil, gen- 
3ll consider 



distribution M MM |23 
entropy [23ll32l [3^361] . coherence visibility: 
eralized purity jS^]. Refs. [H [ll, [H [24 
similar issues on systems comprising more than two sites. 

We mention that Hamiltonian ([T|) can be mapped onto 
a Lipkin-Meshkov-Glick (LMG) model with vanishing 
anisotropy parameter. Some of the results we illustrate in 
the following have been obtained in that context making 
use of numerical or analitical techniques different from 



III. THE SCHRODINGER-LIKE EQUATION 
FOR THE LOW-LYING EIGENSTATES 



Under suitable conditions, the eigenvalue equation for 
Hamiltonian ([1]), H\'i>) = E\'^) is very satisfactorily ap- 
proximated by a Schrodinger-like equation for a fictitious 
particle on a one-dimensional domain. Several slightly 
different versions of this approach have been employed 
the literature [1, [H [13, Hi, IM [H, [iijil. All of them 
are based on the expansion of the eigenstate on the Fock 
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basis of the site occupation numbers 



N 



|vI/)=^c.|7V-i.)i®|i.)2, \v), = ^l=-\Q), (7) 

where |0) is the vacuum state, aj|0) ~ 0. Introducing 
the normalized population imbalance of the Fock state 
relevant to expansion coefficient 



'''' ~ N " ^ " iV 



(8) 



the eigenvalue equation reads 



VlN 



1 + Ziy / 1 — 



1 - Z,, (l + Z 



1 



+ — z^Ct. = i^Ct. (9) 



where the "rescaled energy" E ^ E + U {2N - N^) /4 - 
= i;+ § (iV- 2) - 1] differs from the original eigen- 
value by a unimportant population-dependent constant. 
Now, for large boson populations N, both the square 
roots appearing in in Eq. ^ and the coefficients Cj/±i can 
be expanded in powers of iV~^. Note indeed that the lat- 
ter can be seen as the values that a continuous function 
"0 takes on at the points of the grid mesh in the interval 
[—1, 1] defined by with i/ = 0, 1, 2, • • • ,N. Specifically, 

if we set c^±i = .f^4;{z^±i) = xf^i^iz^ T 2iV-i) 



and retain only the lowest order terms in the expansion of 
the interaction and kinetic term in the eigenvalue equa- 
tion, we obtain the following Schrodinger-like equation 



2 d r d , 



iIj{z) = Eij{z) (10) 



where z € [—1, 1] is the continuum limit of the mesh grid 
and the effective potential is [l^ 



z'^ + e z 



(11) 



Owing to the N~'^ coefficient in front of the derivative 
part in Eq. (|10p . for large populations the low- lying so- 
lutions of such equation will be strongly localized in the 
vicinity of the minima of the potential energy, Eq. (|lip . 
This suggest that the Schrodinger-like equation can 
be further simplified by introducing the (small) devia- 
tion from the minima of Eq. (|lip . Note that the mini- 
mization of Eq. (jlip is equivalent to the determination of 
the lowest-energy stationary solution of the Discrete Sclf- 
Trappin g eq uations arising from the mean-field Hamilto- 
nian ^ |24l | . This is easily verified in the mirror symmet- 
ric case, £ = 0, where, up to the fourth order in the small 
displacement from the minimum point (s), u ^ z — z. 
equation ()10p is equivalent to 



7' 



du^ 

with V'(z) — <p{z — z^) and 



ip{u)=8Lp{u) (12) 
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IV. HARMONIC REGIME 



where 



We note that the leading order in Eq. (fT2|) is almost 
everywhere the second, except at the mean-field critical 
point 7 = — 1 where the first non- vanishing term is the 
fourth. Therefore, for the lowest energy levels, away from 
criticality and for sufficiently large populations, we can 
let = C = and Eq. reduces to the Schrodinger 
equation for a harmonic oscillator, whose discrete spec- 
trum is 



(14) 



1 













-I) 


= ^ ( ^ 


4) 











2VA 



1 



iV|7l\/7^ 



7 > -1 



7 < -1 



(15) 



Actually, for 7 < —1 Eq. describes the small 

oscillations about either of the two equivalent minima 
of the effective potential, z-y = ±y/l + 7~2. Therefore 
a generic solution of Eq. ([T0| is well approximated by 
a superposition of one harmonic solution for each well. 
The symmetry of the problem dictates that these super- 
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positions have a definite parity, 



V2 



(16) 



Also, since in this approximation the two localized func- 
tions (pn{z ± z^) have a vanishing overlap, both ipni^) 
have an energy very close to that in Eq. (fT4)) . In par- 
ticular, the ground-state of the system ipo{z) is well ap- 
proximated by a symmetric superposition "01^(2) of two 
Gaussian functions of the form 



ipa{u) = 



(17) 



with (7^ given in (the lower of) Eq. ((T5|) . 

For 7 > — 1 the effective potential has a single mini- 
mum at 



'7 



0, and the problem can be mapped exactly 
onto a quantum harmonic oscillator, as far as the width 
of the eigenfunctions of the latter problem do not ex- 
ceed the interval in which the harmonic term dominates 
over higher-order terms. The energy spectrum is again 
given by Eq. (jl4p and, in particular, the ground-state is 
iJo{z) = foiz), Eq. ((T7|) . with cr^ given in (the upper of) 
Eq. (HH). 

The most natural physical quantity which is in prin- 
ciple measurable, is given by t he p opulation imbalance 
corresponding to the operator [20| in the original two- 
site problem: 



a\ai — a\a2 
N 



(18) 



Its value and fluctuations are easily evaluated on the 
ground state of the system 



dsz\ijoiz)\^ = 



and kS 



(19) 



(20) 



where (•) denotes the expectation value on the ground 
state of the system. Moreover one can also consider the 
expectation value (•)„ on the nth excited state obtaining: 



(£2) = / dsz^\Mz)\' 



ds z^\Mz)\^ ~ ali2n + 1) + z^^ 



(21) 



The fluctuations of the population imbalance have 
been investigated by several authors, [H, [H, [l^ [23j . 
and analytical expressions for this quantity have been re- 
ported in Refs. [THj, [l^, . In both Gaussian regimes the 
(low lying) energy spectrum of Eq. p2p is well approxi- 
mated by Eq. (Ull). Since = i?„-f[|(Af2_2iV)_l], 
and En depends on as described in table p^ . it is 
easy to calculate the lowest energy gaps in the spectrum 
of H, AEn = En — En^i. As wc mentioned earlier, for 
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FIG. 1: (color online) Low lying energy gaps for A*" = 50 (up- 
per panel) and A'' = 5000 (lower panel). Gray thick lines rep- 
resent the analytic results in Eqs. (|22|) and (I23|l . Specifically, 
the dashed line is Q,^^^ — 1 and the solid line is f^V7 + 1- 
Symbols represent numerical data. Specifically, circles, up- 
ward triangles, squares and downward triangles are Ai5„ with 
n = 1, 2, 3, 4, respectively. 



71 <C 7 <C 72 each energy level is (almost) twofold de- 
generate, so that we have 



AS. 



2n+l 



0, 



7i < 7 < 72 



+ 1, 73 < 7 < 74 



and 



AEo 



1, 



71 < 7 < 72 
73 < 7 < 74 



(22) 



(23) 



In Fig. [T]we compare Eqs. and ([25)1 with numerical 
data for A£'„. The plots clearly show that the analytic 
approximations become more accurate with with increas- 
ing boson population. 



V. LIMITS OF VALIDITY AND REGIMES 

In the present section we describe the limits of valid- 
ity of the harmonic solutions described above and the 
different regimes taking over at different values of the 
effective parameter 7. A crucial assumption in passing 
from Eq. ^ to Eq. PH)) is that the coefHcients can be 
regarded as continuous functions of the population im- 
balance of the Fock state defined in Eq. ([5]). This may 
be problematic for eigenstates in the mid-spectrum of 
Hamiltonian ([1]), but is very reasonable for states close 
to the extrema of the same spectrum. However, the con- 
tinuous approach can prove ill posed even for the ground 
state of the system. It is clear that, for the continuous 
limit of Eq. (|10p to be well defined, the Gaussin width 
must be much larger than distance between two neigh- 
boring points in the discrete Eq. ©, namely 2N~^. This 
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means 7 » 71 ^ 
for 7 > -1 [ii]. 



-Vn for 7 < -1 and I7I < 74 iV^ 



We observe that the same estimates can be obtained 
by imposing that the expansion parameter in the strong 
couphng regimes be actually small, as discussed in some 
detail in Sec. [A] This ensures that for 7 ^ 71 and 7 ^ 
74 the first-order perturbativc description of the ground 
state applies. 

The parameter range 71 <C 7 74 where the continu- 
ous approach applies can be divided into three regimes. 
For 71 <C 7 <C 72 the ground state of the system is well 
approximated by a symmetric superposition of two Gaus- 
sians. On the other hand, for 73 <C 7 <C 74 the ground 
state is well approximated by a single Gaussian. 

These two regimes are separated by an intermediate 
critical regime 72 < 7 < 73, surrounding the mean- 
field critical point 7c = —1, where the higher order 
terms in Eq. (|T2|) cannot be neglected and the prob- 
lem is not Gaussian any more. The boundaries of this 
regime can be obtained from a comparison between the 
leading and next-to-leading term in the pcrturbative ex- 
pansion of the effective potential. A reasonable esti- 
mate for the lower bound 72 < 7c can be obtained ei- 
ther by imposing yl cr^ |B|cri^ or by requiring that the 
distance between the two equivalent harmonic minima, 
2yJ\ — 7~^, is much larger than the width of the rele- 
vant Gaussian ground states . In both cases one finds 
7c — 72 ~ N~'^/^ . Similarly, the upper bound 73 > 7c is 
obtained by requiring ^ct^ ^ Ccri, which results once 
again into 73 — 7c ^ N~^^^ (Tol . |45| . As discussed in 
Sec. IVIl the energy gaps in the critical regime depend on 
the system population, AE ^ N~^^^, at variance with 
what happens in the Gaussian regimes. This provides a 
iV— dependent cutoff to the vanishing of AE at 7 = 7c 
predicted in Eqs. ([22]) and ((23)) . 

In summary, one can recognize five different regimes 
depending on the value of the effective parameter 7. For 
7 <C 71 and 7 ^ 74 the continuous approximation of Scc- 
tion lllll does not apply, but the system can be analyzed by 
resorting to first-order perturbativc theory, as discussed 
in Appendix [X] For 71 <C 7 ^ 72 the continuous ap- 
proximation applies, and the ground-state of the system 
is a symmetric superposition of two Gaussians centered 
at the equivalent harmonic minima of the effective po- 
tential. For 72 ^ 7 <IC 73 the system is critical, in a 
sense that will be clarified shortly. Note that in the limit 
of large populations the critical regime reduces to the 
mean-field critical point. For 73 ^ 7 ^ 74 the ground- 
state is a Gaussian centered at the unique minimum of 
the effective potential, = 0. Note that in estimating 
the location of the boundaries 71-4 we only considered 
the dependence on the (large) bosonic population, and 
neglected other numerical factors. 



VI. CRITICAL REGIME 

According to our previous analysis, as the effective 
parameter 7 approaches the mean-field critical value 
7c = —1, the system exhibits the hallmarks of a phase 
transition, where the minimum Zj of the effective po- 
tential, table ([13]), plays the role of the order parame- 
ter. Note that this quantity is in principle measurable as 
the absolute value of the population imbalance operator 
Zj = ( 1 2 1) . The fluctuations of this order parameter 

(|zp)-(|z|)2=a^ (24) 

coincide with the Gaussian width of Eq. (fT5|) . 



The behavior for z^ and a:^ described in table ([T^ and 
Eq. (j24p , respectively, holds true in the mirror-symmetric 
case, V = e = 0. It proves useful to study behaviour of 
the order parameter in the presence of a vanishing energy 
offset between the sites of the dimer, \v\ ^ 1. According 
to our discussions about the Gaussian regimes, the value 
of the order parameter will coincide with the absolute 
minimum of the effective potential which, owing to the 
presence of a non- vanishing e, is now non degenerate for 
any value of 7. The minimization of Eq. (jlip then gives 



(25) 



where only linear terms in e have been taken into account. 
Eq. ([25l) allows us to introduce a sort of susceptibility 



Xz 



de 



{m)e = 




(26) 



Equations (gll), (EH) provide some critical 

exponents for the system. As 7 approaches its critical 
value 7c = —1, the order parameter vanishes (from be- 
low) as (|z|) ^ (7c — j)"' , its fluctuations diverge as 
— (|z|)^ ~ l7^7c|^"" and its susceptibility diverges 
as Xz ^ I7 ~ 7c I where 



I 

2' 



1 



(27) 



The above discussion rigorously applies in the limit of 
inflnite population, which plays the role of an effective 
thermodynamic limit. Indeed, as we observe in Section 
IVl only in this limit does the critical regime shrink to a 
single point corresponding to the critical value. For any 
finite populations there exists a finite region surrounding 
the critical point where the Gaussian results do not apply. 
In this region the quartic term dominates and Eq. (jlOp 
can be approximated as 



dz^ 



Cz" 



(28) 



where C is listed in the upper row of table ([T3| . It is easy 
to prove that in this critical regime the eigcnfunctions of 
Eq. dUl) have form 



V'„(z) = N'/'^MzN'/^), 



(29) 
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FIG. 2: (color online) Dependence on the population of the 
first few energy gaps AE„ at 7 = — 1 (from bottom to top, 
n = 1, 2 • • • , 5). The data points are numerically determined 
values, the solid lines are just guides to the eye. 



where 4>niC) is the solution of the population- 
independent problem 



16 



4>niC) = £n4>niC) 



(30) 



and the eigenvalue corresponding to ipniz) is £„ = 
N^^^Sn- This result allows us to estimate the behaviour 
of the (low-lying) energy gaps in the spectrum of Hamil- 
tonian ([T]) in the critical regime. We get AEn ^ N^'^^^, 
where we used the definition of £ in table (|13p and the 
fact that the spectrum of problem (pO| does not depend 
on N. This estimate, whose correctness demonstrated in 
Fig. [21 provides an alternate route for the determination 
of the boundaries of the critical region. Indeed, these can 
be assessed by requiring that the prediction in the Gaus- 
sian regimes, Eqs. ([^^ and are of the same order 
as that in the critical regime, i.e. This gives once 

again 7c — 72 iV^2/3 g^^^^ „ ^ ^-2/3^ 

As to the fluctuations of the order parameter, partic- 
ularizing Eq. ([^ to the ground-state of the system, 
ri = 0, we get 

jdce\Mc)\' 



N' 



-2/3, 



JdClMCW 



(31) 



Within a standard scaling approach it is reasonable to 
assume that there exists a correlation length, dictating 
the effective size of the system, which diverges at the 
critical point 



Af-y 



17 -7c 



(32) 



Any physical quantity should depend on the boson pop- 
ulation N only through the ratio N/J\f^, so that it should 
be 



{\z\')-{\z\r^N^f 



N 



(33) 
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FIG. 3: (color online) Scaling of the fluctuations of the (ab- 
solute value of) the population imbalance operator. 



The comparison with Ec^s. (|15p . ([M)) and (pij) then allows 
us to conclude that 



(34) 



3' ^ 2 



The correctness of the scaling assumption in Eqs. (|33p - 
(pl)) is demonstrated in Fig. [3l where a nice data collapse 
is observed for a wide range of boson populations. We 
mention that - in the context of the LMG model - the 
critical exponents in (j34p were determined numerically in 
Ref. jssl l and making use of an analytic approach different 
from ours in Refs. [39l44l| . 

We conclude by remarking that the rescaled energies 
£n and the expectation values for the rescaled problem 
(pO)) can be estimated by applying the Bohr-Somnicrfeld 
quantization rule. This gives 



2£i'- r ^ 1 

i/^V^"" 16'^^^^^"+2^' 



2£ 



(35) 



which provides quite satisfactory estimates of the energy 
levels of the system, 



i?„ = J7CiiV-V3(n + 1)2/3 



(36) 



where Ci is a suitable constant. Within the same ap- 
proximation the mean square displacement of the eigen- 
function (/'n(C) can be evaluated as well. In particular, 
for the simclassical orbit of energy £„, one gets 



(37) 



where C2 is a suitable constant. 
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VII. QIT INDICATORS 

As we mention above, quantities borrowed from quan- 
tum information theory have been employed to investi- 
gate the the incipient quantum phase transition between 
the locahzed and delocahzed ground-state of the system. 
Owing to the particularly simple structure of the system, 
some of these quantities turn out to be entirely equivalent 
to usual indicators. This is the case of the generalized 



purity of Ref . [37[ , which basically coincides with the co- 
herence visibility considered in Refs. HEHEHIIj i-e- the 
off-diagonal term of the one-body density matrix. Simi- 
larly, the Fisher information of Ref. is equivalent to 
the population imbalance fluctuations considered in sev- 
eral earlier works [HI, [l^ [l^ [2^ as well as in the present 
paper. 

The fidelity and entanglement entropy considered in 
Refs. 31 1 and [l^, respectively, appear to be less trivial 
quantities [49j . The former is nothing but the overlap 
between two ground-states of the system corresponding 
to slightly different values of the parameter driving the 
transition (soj . Since this quantity is almost 

everywhere extremely close to unity, a more convenient 
indicator is provided by the so-called fidelity susceptibility 

M 



lim 



1 



(*7l*7+0 



2dK^ 



1=0 



(38) 




-1.5 



-0.5 

7 



0.5 



FIG. 4: (color online) The analytic expression obtained for 
the fidelity susceptibility in the Gaussian regimes is compared 
with the numerical results for the same quantity. The circles, 
upward triangles, squares, downward triangles and diamonds 
are obtained from exact diagonalization of Eq. ([T]) for sys- 
tems containing 500, 1000, 2000, 4000 and 8000 bosons, re- 
spectively. In all cases we used Eq. (|38|) with k — 10~^. The 
thick gray lines represent Eq. (|40p for the same boson popu- 
lations. Notice how Xf(7) depends on the boson population 
only for 7 < 7c. The inset contains a magnification of the 
region enclosed by the dashed black frame, where Eq. (|40p 
does not apply. The thin solid lines are guides to the eye. 



In a bipartite system, the entropy of entanglement of a 
pure state \^) is given by — Tr(pj- log pj), where pj is the 
reduced density matrix of subsystem j, obtained from the 
full density matrix after tracing over the states of 

the other subsystem [52] ■ In the case under investigation 
the two subsystem are of course the two sites composing 
the Bose-Hubbard dimer and, owing to number conser- 
vation, the entropy of entanglement of an eigenstate of 
Hamiltonian ([IJ is 



1 



log2(A^+l) 



JV 

■El 



'logzlc^l^ 



(39) 



where the Ci/'s are the expansion coefficients in Eq. (O 
and we introduced a normalization factor. This quan- 
tity has been considered in a few earlier works on the 
Bose-Hubbard dimer. In the presence of repulsive inter- 
actions it was used to investigate the precursor of the 
Mott insulator-superfluid quantum phase transition oc- 
curring on infinite lattices [ssj . In attractive systems 
it was employed for the characterization of Schrddinger- 
cat-like states in the presence of a small energy offset v 
between the two sites of the system [s^ , and in the study 
of the quantum counterpart of the mean-field transition 
for small finite boson populations [23j . 

These two indicators have been considered for the 
LMG model [l^, HJ] which, as we mention, can be 
mapped onto Hamiltonian ([IJ. However, the natural de- 
composition of the system necessary for the calculation 



of the entanglement entropy is quite different in the spin 
and boson context, and the results are not readily com- 
parable. As to the ground-state fidelity, a comparison 
will be made with the results in Ref. j43 |. 

Based on the results sketched in the previous sections 
and in the appendices, we are able to give an analytic 
description of the behaviour of the fidelity susceptibility 
and entropy of entanglement, Eqs. (|55)) and ([M)) . 

Let us start with the former quantity. As we discuss 
in Section Hill in the Gaussian regimes the ground-state 
of the system is strictly related to the ground state of 
a harmonic oscillator. Specifically, for 71 ^ 7 ^ 72 
the ground state is well approximated by a symmet- 
ric superposition , ip^lz), of two Gaussian functions 



of square width cr. 



centered at 



= ±-\/l — 7^2, as described by Eqs. (fTH]) and ()17p . 
Likewise, for 73/I7 <ti 74 the ground state is well de- 
scribed by a single Gaussian ipo{z) of square width cr^ = 

{Ny/1 + 7) ^ centered at z-y = 0, as described by Eq. 
pri) . Straightforward calculations result in 



JV 



Xf(7)= 



+ 7i«7«72 ^^^^ 

^ 73 ^ 7 <C 74 



64(7-1-1)^ 



whose correctness is demonstrated in Fig. 2] It is 
worth noticing that Eq. PO)) predicts a dependence on 
the boson population only in the lower Gaussian regime 
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7i 7 <C 72, whereas for 73 ^ 7 <C 74 the fidehty 
susceptibihty turns out to be A^— independent. This can 
be explained by the fact that in the latter regime the 
small variation k in the effective parameter affects only 
the width of the single Gaussian representing the ground- 
state of the system. Conversely, for 71 <C 7 ^ 72, the 
change affects both the centers and the widths of the 
two Gaussians comprising the ground state of the sys- 
tem. Also notice that Eq. pO|) predicts a divergence in 
the fidelity susceptibility at the mean- field critical point, 
7 = 7c. However, for any finite N , there is a small region 
surrounding such point where the Gaussian approxima- 
tions do not apply. This is demonstrated in the inset 
of Fig. m which clearly shows that for finite populatons 
Xf(7) features a finite peak in the vicinity of 7c. The 
larger the population, the sharper and the closer to the 
critical point is the peak. 

In order to study the behavior of xf(7) in the critical 
region it proves useful to exploit Eq. ([25)) . As it is dis- 
cussed in Ref. [s^, the limit in Eq. ([55)1 can be carried 
out exactly through a perturbative expansion involving 
the entire spectrum of the problem at a given value of 
the driving parameter 7. Considering the harmonic term 
in Eq. (|28p as a perturbation over the quartic term, we 
get 



Xf(7c) 



16 



E 



/dCC'0o(O0n(C) 



£0 ~ 



(41) 



where the eigenfunctions (t>n{0 ^-nd eigenvalues are 
dfined by (jSO]) and (|30l) . Rigorously speaking the 4>n{C) 
provides a good approximation for an eigenstate of the 
original problem, Eq. ([T]), only for sufficiently small val- 
ues of the quantum number n. However, only the first 
few terms in the sum appearing in Eq. (|4ip are expected 
to provide a significant contribution, owing to the lo- 
calization of the ground state 4'aiO and the energy gap 
appearing in the denominator. 

The N'^/^ behavior predicted in Eq. (|4ip and con- 
firmed by the numerical results in Fig. [5] can be rec- 
ognized as the superextensive divergence of the fidelity 
susceptibility which, according to Ref. [5l|, is the hall- 
mark of a quantum phase transition. The same exponent 
was estimated numerically in Ref. [l^]. Also, our Eq. PO)) 
agrees with a similar result in the same paper, although 
we obtain a different behavior of the sub-leading terms. 

Let us now turn to the entropy of entanglement. Re- 
calling that the continuous version of Eq. ([39)) is 



-^/_^dz|^(z)|^log,(||^(z)n, (42) 



it is easy to check that in the Gaussian regimes 

1 , C{j) 



Seil) 



2 21og2iV 



(43) 
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FIG. 5: Scaling of the fidelity susceptibility at the mean- 
field critical value. The data points are obtained from exact 
diagonalization of system containing A'' — 500 ■ 2" bosons, 
with n = 0, • ■ • ,5. The thin solid line corresponds to 7.1 ■ 
10-3 Ar4/3. 



where 



log2 



C(7) = ( \^h\V^i 



71 < 7 < 72 



1082(2;^) 73«7«74 



(44) 



Note that, despite Eq. (|39)) forbids that S'e(7) exceed 
unity, the quantity in Eq. ()43p diverges in the vicinity of 
the mean-field critical point. However, this is not alarm- 
ing, since at any finite boson population there is a small 
interval around such point where the ground-state is not 
Gaussian and Eq. (|^5)) does not apply. It could be easily 
proven that the requirement that Eq. ([43| be less than 
unity provides an estimate for the boundaries of the non- 
Gaussian regime that is equivalent to the one derived in 
Section |V| An interesting feature revealed by Eqs. (|^5)) - 
(|44)) is that in the — J- 00 limit the entanglement en- 
tropy tends to a constant, independent of the interaction 
to hopping ratio, £'0(7) = ^ V7 7^ 7c. Note indeed that 
quantity in Eq. (j^^ does not depend on TV, so that the 
second term in Eq. (|43p vanishes with increasing popula- 
tion. 

The behavior of the entropy of entanglement in the 
critical regime can be obtained by plugging iIj{z) = 
N^/^M^^^^z) into Eq. (|42]), where MC) 

is once again 

the ground state of Eq. (|5(I)) . Straightforward calcula- 
tions give 



^c(7) = o 



1 



dC|0o(C)Plog2 [2|0o(C)l'] (45) 



3 log2 N 



Since the integral in the second term of Eq. (j^ does 
not depend on N, we get that limw-s-oo 'S'e(7c) = §■ Thus 
the mean-field critical point is signalled by the entropy 
of entanglement as a peak which, in the large popula- 
tion limit, turns into a discontinuity. Figure [5] shows 
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FIG. 6: (color online) The analytic expression for the entan- 
glement entropy, Eq. (|43p . is compared with numerical results 
at different boson populations. Numerical data (symbols) and 
analytic results (solid lines) corresponding to the same boson 
population have the same color. The black dashed horizontal 
line is the limit of the Gaussian result, Eq. ((43]). The black 
star is the thermodynamic limit of the entanglement entropy 
at the mean-field critical point, Eq. (|45|l . The vertical dashed 
line signals the mean-field critical point. 



a comparison between Eqs. (|43p - (|44p and numerical re- 
sults obtained via exact diagonalization. The agreement 
is quite satisfactory, and improves with increasing boson 
population. Note that even at relatively large popula- 
tions the entropy of entanglement is still rather different 
from its limit. This is because the finite size corrections 
decrease logarithmically with N. The different behavior 
at criticality is demonstrated in Fig [7l which shows the 
correctness of Eq. (|45p . 



VIII. THERMAL FLUCTUATIONS 

We now turn to the thermal fluctuations of the order 
parameter z 



(46) 



where {■)t and (•)„ denotes the thermal average and the 
quantum expectation value on the state of energy En, 
respectivley. 

Evaluating i?„, (z^)„ and (z)„ using standard proper- 
ties of quantum hamonic oscillators (see Section ITV]) we 
get 



(47) 



FIG. 7: Scaling of the entanglement entropy at the mean- 
field critical value. The data points are obtained from exact 
diagonalization of system containing A*" = 50 ■ 5" bosons, with 
71 = 0, • • • ,5. The thin solid line corresponds to 3.4 logj A'^ — 
0.11. 



where the thermal fluctuations are 



coth 



N\f\y/j^~l 



7i < 7 < 72 



73 < 7 < 74 



(48) 



and Z7 is given in Table (fT3|) . Therefore for U-qT < AE2 
(i.e. the lowest energy gap of the system) a^ ,^ is es- 
sentially given by the fluctations of the quantum ground 
states ([20]) . whereas for large temperatures the thermal 
fluctuation dominates and we have 



nWM^' 7i«7«72 



2fcBT 
OAr(7-|-l) 



(49) 



73 < 7 < 74 



Eq. highlights the typical linear dependence on 

temperature of the thermal fluctuations in harmonic sys- 
tems. We remark that both thermal and quantum fluc- 
tuations decrease with the total population as 1 /TV. This 
is demonstrated by the data collapse in Fig. |S1 

Also in the thermal regime we expect the harmonic ap- 
proximation to apply only when the higher order therms 
in Eq. ((T^ arc negligible. This means that for 71 < 



7 < 72, o'^ ^ should be smaller than A/B while for 
73 < 7 < 74, (t|^^ should be smaller than A/C. In both 
cases for 7 w —1 we get kT < N{j + 1)^. Note indeed 
that, as evident from Fig. [8l the higher the temperature, 
the larger must the population be for the harmonic ap- 
proximation (dash-dotted line) to apply. 

Let us now consider the regime 72 < 7 < 73 where 
the quartic term dominates. The values of £„, {z)t and 
{z^)n can be estimated within the semiclassical Bohr- 



Sommerfeld approximation discussed in Sec. I VII Plug- 
ging Eqs. (1251) and into Eq. (gH) we get, for 
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FIG. 8: (color online) Scaling of imbalance fluctuations as 
a function of temperature for 7 = —0.2. Black, red, blue 
and magenta solid line correspond to 100, 200, 1000 and 2000 
bosons, respectively. Dashed and dash-dotted line is the the- 
oretical predictions (|15|l and (|49p . respectively. 



7 = 7c = -1, 



(50) 

It is easy to check that also in this critical regime ct^^t 
reduces to the purely quantum contribution when fceT' is 
smaller than the energy gap between the lowest eigenpair, 
Cifl/N^^'^ . At higher temperatures we have 



2n 



3 dn 



iVi/e" 



7173' 



' dn 



kBT 



(51) 



which is the typical behaviour of classical thermal fluc- 
tuations when a pure quartic potential is present. We 
remark that within the critical region the fluctuations 
of the population imbalance due quantum effects scale 
as N~'^/^ with increasing population, while the thermal 
fluctuations scale as N~-^/^. This is shown in the upper 
and lower panel of Fig. |9l respectively. 

The different quantum and classical regimes around 
7 = — 1 can be summarized as in Figure 1101 Note that 
the boundaries between different regions, signalled by the 
coloured lines, are actually crossovers between different 
behviours. The only true critical point is 7 = —1, in the 
limit of infinite N. 

In regions I and III cr^ ^ coincides with the result ap- 
plying at zero-temperature for purely quartic and har- 
monic potentials, respectively. Region II is characterized 
by thermal fluctuations in a harmonic potential, while 
in region IV we have a classical behaviour in an anhar- 
monic potential. In regions I and II cr^ j, decays as 7V~^, 
while in region III quantum critical effects produce fluc- 
tuations decaying as N~'^/^. Finally, in region IV fluc- 
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N = 1 00 
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kBT/n 

FIG. 9: (color online) Scaling of imbalance fluctuations as a 
function of temperature for 7 = 7c = — 1 and different bo- 
son populations. The same data are multiplied by a different 
function of the boson population in the two panels, to high- 
light the different scaling laws in the quantum and thermal 
regime. The thick gray line in panel B is the typical behaviour 
of the classical thermal fluctuations for a pure quartic poten- 
tial, Eq. ((51]). 
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FIG. 10: (color online) In this figure we represent the different 
thermal and quantum regimes in the region aroud the critical 
point 7 = — 1 the data are obtained for TV = 100. In regions I- 
IV represent quantum harmonic, thermal harmonic, quantum 
quartic and thermal non-harmonic behaviours respectively. 



tuations should decay as ^^'^ for not too large tem- 
peratures, and they should be independent of N for very 
large kBT/Q. 



IX. CONCLUSIONS AND PERSPECTIVES 

In this paper we addressed the quantum counterpart 
of the transition from a localized to a delocalized ground- 
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state occurring in an interacting bosonic gas trapped in a 
double well potential, i.e. a bosonic Josephson junction. 
Wc identified and discussed five regimes depending on 
the value of the effective interaction among bosons. For 
sufficiently small effective interactions we analyzed the 
low-lying spectrum of the system Hamiltonian (a two- 
mode Bosc-Hubbard model) by recasting the problem in 
terms of a Schrodinger-likc equation for a single particle 
in a one-dimensional domain. In particular we studied 
how the results of the analysis on the quantum system 
reduce to those of the semiclassical treatment based on 
the Discrete Self- Trapping equations as the bosonic popu- 
lation is increased. We devoted a special attention to the 
critical regime which, in the theromodynamic limit of in- 
finite bosonic population, collapses to a bifurcation point 
in the semiclassical picture. We extended our analysis to 
quantities borrowed from Quantum Information Theory, 
namely the entropy of entanglement and the ground-state 
fidelity, and to finite temperature effects. In the latter re- 
spect, we evidenced that, in the critical regime, the quan- 
tum and thermal fluctuations of the population imbal- 
ance of the two wells exhibit different scaling behaviors, 
which may be amenable to quantitative measurement. 
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Appendix A: Strong-coupling regimes 

As we mention in Section |Vl when interaction domi- 
nates over kinetic energy the continuous approach of Sec- 
tion mil does not apply. However, in this situation an 
approximation of the ground-state of the system can be 
given by resorting to perturbative theory. For large re- 
pulsive interactions, first order perturbative theory gives 

|vI/) = |vI/o)+ A (|vl/+) + !*_)) (Al) 
where, assuming without loss of generality that N is even. 



from the "single Gaussian" regime, 73 ^ 7 <C 74, to 
the present "strong coupling" regime, 7 ^ 74, there is a 
change in the behaviour of the fluctuations of the popu- 
lation imbalance operator, Eq. (jlSp . on the ground-state 
of the system, as discussed in Ref. [2^. Indeed in the 
former regime this is clearly (z^) — (z)^ = ct^. where 
(T^ ^ A^"^ is given in Eq. (|15p . Conversely, in the strong 



(A4) 



couphng regime of perturbative state l|Aip . we get 



8A2 



7V2(l + 2A2) 



2^ 



where the last approximate equality applies in the limit 
of small perturbations and large populations. Since for 
7 ~ 74 the fluctuations in Eqs. (|15p and (jA4p arc of 
the same order of magnitude, there is no intermediate 
regime between the "single Gaussian" and the repulsive 
"strong-coupling" regimes 

Similar arguments apply in the attractive "strong- 
coupling" regime. At first order in perturbation theory 
the ground state is 



|*) = |*o)-^(|*+) + |*-)) (A5) 



where Ivfo) 



1 

V2 



|l,iV- 



(1^,0) 
1) and A = 



|0,iV)), 1*+) = |7V-1,1), 



The reqircmcnt that 



2U N-1 ■ 

A <C 1 identifies the attractive "strong-coupling" regime 
7 <C 71 ^ —Vn. Straightforward calculations give 



1 + A 



2 {N-2f 



Ar2 



1+A2 



r 



(A6) 



We observe that, unlike the repulsive case, there is no 
cross-over of the fluctuations of operator in passing 
from the strong-coupling regime 7 ^ 71 to the "double 
Gaussian" regime described in Section [V] Note indeed 
that Eq. (|A6p coincides with Eq. (|20|) at large negative 
7. The same results apply when one considers a single 
peak. 



and 



N N 
2'^ 



N N 

— ± 1, — T 1 
2 ' 2 ^ 



(A2) 



(A3) 



Equation (|Aip applies when the factor in front of the 
first order perturbative term is small, A ^ 1. It is easy 
to check that, for this to be true, it should be 7 74, 
where 74 ~ N"^ is the upper bound of the "single Gaus- 
sian" regime discussed in Section [V] Note that in passing 



I*) = |iV,0) + A|Ar- 1,1) 



which gives 



1 + A^ 



1 + A2 



1 + 
1 + A2 



(A7) 



(A8) 



But this is exactly the same behavior as cr^ in Eq. (jlSp 
for large populations and large (negative) 7. 
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